﻿using System;
using System.Numerics;
using System.Threading.Tasks;

namespace Mrkz.ImageManipulation.FrequencyFiltering
{
    /// <summary>
    ///
    /// </summary>
    internal static class FFTEngineSupportFunctions
    {

        internal static int calc_pow(int h)
        {
            int p = 1;
            while (p < h)
            {
                p *= 2;
            }
            return p;
        }

        public static TImage<ComplexPixel> ParallelFFT2D(TImage<ComplexPixel> InputData, int Direction)
        {
            int Height = calc_pow(InputData.Height);
            int Width = calc_pow(InputData.Width);

            Complex[,] out_R = new Complex[Height, Width];
            Complex[,] out_G = new Complex[Height, Width];
            Complex[,] out_B = new Complex[Height, Width];

            //per ogni colonna            
            Parallel.For(0, Width, (x) =>
            {
                Complex[] _R = new Complex[Height];
                Complex[] _G = new Complex[Height];
                Complex[] _B = new Complex[Height];

                //copio tutta la colonna
                for (int y = 0; y < Height; y++)
                {
                    ComplexPixel p = InputData.GetPixelWithPad(y, x);

                    _R[y] = p.CRed;
                    _G[y] = p.CGreen;
                    _B[y] = p.CBlue;
                }

                _R = NormalizedFFT1D(_R, Direction);
                _G = NormalizedFFT1D(_G, Direction);
                _B = NormalizedFFT1D(_B, Direction);

                //copio la colonna nell'array output
                for (int y = 0; y < Height; y++)
                {
                    out_R[y, x] = _R[y];
                    out_G[y, x] = _G[y];
                    out_B[y, x] = _B[y];
                }
            });


            //per ogni riga
            Parallel.For(0, Height, (y) =>
            {
                Complex[] _R = new Complex[Width];
                Complex[] _G = new Complex[Width];
                Complex[] _B = new Complex[Width];

                //copio tutta la riga
                for (int x = 0; x < Width; x++)
                {
                    _R[x] = out_R[y, x];
                    _G[x] = out_G[y, x];
                    _B[x] = out_B[y, x];
                }

                _R = NormalizedFFT1D(_R, Direction);
                _G = NormalizedFFT1D(_G, Direction);
                _B = NormalizedFFT1D(_B, Direction);

                //copio la riga nell'array output
                for (int x = 0; x < Width; x++)
                {
                    out_R[y, x] = _R[x];
                    out_G[y, x] = _G[x];
                    out_B[y, x] = _B[x];
                }
            });

            return ImageFactory<ComplexPixel>.FromComplexPixelArray(Height, Width, InputData.CropArea, out_R, out_G, out_B);
        }

        public static TImage<ComplexPixel> FFT2D(TImage<ComplexPixel> InputData, int Direction)
        {
            int Height = calc_pow(InputData.Height);
            int Width = calc_pow(InputData.Width);

            Complex[,] out_R = new Complex[Height, Width];
            Complex[,] out_G = new Complex[Height, Width];
            Complex[,] out_B = new Complex[Height, Width];

            //per ogni colonna            
            for (int x = 0; x < Width; x++)
            {
                Complex[] _R = new Complex[Height];
                Complex[] _G = new Complex[Height];
                Complex[] _B = new Complex[Height];

#if TIMING
                TimeSpan t1,t2,t3;
                System.Diagnostics.Stopwatch s = new System.Diagnostics.Stopwatch();
                s.Start();
#endif
                //copio tutta la colonna
                for (int y = 0; y < Height; y++)
                {
                    ComplexPixel p = InputData.GetPixelWithPad(y, x);

                    _R[y] = p.CRed;
                    _G[y] = p.CGreen;
                    _B[y] = p.CBlue;
                }

#if TIMING
                t1=s.Elapsed;
#endif
                _R = NormalizedFFT1D(_R, Direction);
                _G = NormalizedFFT1D(_G, Direction);
                _B = NormalizedFFT1D(_B, Direction);
#if TIMING
                t2 = s.Elapsed;
#endif
                //copio la colonna nell'array output
                for (int y = 0; y < Height; y++)
                {
                    out_R[y, x] = _R[y];
                    out_G[y, x] = _G[y];
                    out_B[y, x] = _B[y];
                }
#if TIMING
                t3 = s.Elapsed;
                s.Stop();
                Console.WriteLine("{0}ms {1}ms {2}ms",t1.TotalMilliseconds,t2.TotalMilliseconds,t3.TotalMilliseconds);
#endif
            }


            //per ogni riga
            for (int y = 0; y < Height; y++)
            {
                Complex[] _R = new Complex[Width];
                Complex[] _G = new Complex[Width];
                Complex[] _B = new Complex[Width];

                //copio tutta la riga
                for (int x = 0; x < Width; x++)
                {
                    _R[x] = out_R[y, x];
                    _G[x] = out_G[y, x];
                    _B[x] = out_B[y, x];
                }

                _R = NormalizedFFT1D(_R, Direction);
                _G = NormalizedFFT1D(_G, Direction);
                _B = NormalizedFFT1D(_B, Direction);

                //copio la riga nell'array output
                for (int x = 0; x < Width; x++)
                {
                    out_R[y, x] = _R[x];
                    out_G[y, x] = _G[x];
                    out_B[y, x] = _B[x];
                }
            }

            return ImageFactory<ComplexPixel>.FromComplexPixelArray(Height, Width, InputData.CropArea, out_R, out_G, out_B);
        }

        /// <summary>
        /// Implementazione della FFT secondo cooley-tukey su array di numeri complessi
        /// </summary>
        /// <param name="InputArray">input data</param>
        /// <param name="Direction">+1 forward, -1 reverse</param>
        /// <returns>result array</returns>
        public static Complex[] NormalizedFFT1D(Complex[] inputarray, int dir)
        {
            Complex[] X = FFT1D(inputarray, dir);
            double amplitude = Math.Sqrt(X.Length);

            for (int i = 0; i < X.Length; i++)
                X[i] /= amplitude;

            return X;
        }

        /// <summary>
        /// Implementazione non normalizzata della FFT 
        /// </summary>
        /// <param name="InputArray"></param>
        /// <param name="Direction"></param>
        /// <returns></returns>       
        public static Complex[] FFT1D(Complex[] InputArray, int Direction)
        {
            if (InputArray.Length == 1)
                return InputArray;

            int N = InputArray.Length;
            Complex[] OutputArray = new Complex[N];
            Complex[] OddData, EvenData;

            EvenData = new Complex[N / 2];
            OddData = new Complex[N / 2];
            for (int k = 0; k < N / 2; k++)
            {
                EvenData[k] = InputArray[2 * k];
                OddData[k] = InputArray[2 * k + 1];
            }

            OddData = FFT1D(OddData, Direction);
            EvenData = FFT1D(EvenData, Direction);

            double ph = Direction * -2 * Math.PI / N;

            for (int k = 0; k < N / 2; k++)
                OddData[k] *= Complex.FromPolarCoordinates(1, ph * k);

            for (int k = 0; k < N / 2; k++)
            {
                OutputArray[k] = EvenData[k] + OddData[k];
                OutputArray[k + N / 2] = EvenData[k] - OddData[k];
            }

            return OutputArray;
        }



    }

}
